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Abstract:  An  extended  Kalman  filter  was  developed  to  values  in  a  log-transformed  metric.  At  St.  John  River  at 

automate  the  real-time  projection  of  ice-affected  Dickey,  Mainejogarithms  of  projected  streamflow  values 

streamflow,  based  on  routine  measurements  of  stage  were  within  8%  of  the  logarithms  of  published  values 

and  air  temperature  and  the  relation  between  stage  87.2%  ofthe  time  and  within  15%  of  published  values 

and  flow  during  open-water  conditions.  The  form  96.6%  of  the  time  during  periods  of  ice  effects.  At  Platte 

accommodates  three  dynamic  modes  of  ice  effects:  River  at  North  Bend,  Nebraska,  logarithms  of  projecfed 

sudden  formation-ablation,  stable  ice  conditions,  and  streamflow  values  were  within  8%  of  the  logarithms  of 

final  elimination.  The  filter  was  applied  to  historical  published  daily  values  90.7%  ofthe  time  and  within 

data  from  two  long-term  streamflow-gaging  stations.  15%,  97.7%  of  the  time  during  ice-affected  condi- 

They  were  stable  and  parameters  converged  for  bofh  tions.  This  extended  Kalman  filter  allows  estimation  of 

stations,  producing  estimates  that  were  highly  corre-  ice-affected  streamflow  at  other  gaging  stations  by  ad- 

lated  with  and  linearly  related  to  published  streamflow  justing  filter  parameters  to  site-specific  conditions. 


Coven  Ice  jam  in  1991  atSf.  John  River  at  Dickey,  Maine,  gage. 
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Projecting  Ice- Affected  Streamflow 
by  Extended  Kalman  Filtering 

DAVID  J.  HOLTSCHLAG,  CHARLES  T.  PARKER, 
AND  MOHINDER  S.  GREWAL 


INTRODUCTION 

The  U.S.  Army  Corps  of  Engineers  is  responsible  for  mitigating  flood  damage  and  maintaining 
navigable  conditions  in  rivers  throughout  the  U.S.  Decisions  affecting  the  day-to-day  operation  of 
Corps  dams  are  based  on  real-time  streamflow  data  obtained  by  telemetry  from  streamflow-gaging 
stations  typically  operated  by  the  U.S.  Geological  Survey  (USGS).  Ice  effects — variable  blockage  of 
the  channel  by  ice  and  increased  flow  resistance — reduce  the  accuracy  of  real-time  streamflow  data. 
To  improve  the  controllability  of  flow  during  ice-affected  periods,  CRREL  supported  this  effort  to 
improve  real-time  estimates  of  ice-affected  streamflow. 

USGS  operates  a  network  of  more  than  7200  continuous-record  streamflow-gaging  stations 
nationwide  (Condes  de  la  Torre  1994).  Direct  measurements  of  flow  and  stage  (water-surface  elevation 
in  the  waterway)  are  routinely  obtained  about  every  6  weeks.  The  average  relationship  between  stage 
and  flow  during  open-water  (no  ice  cover)  conditions — the  stage-streamflow  rating — monotoni- 
cally  increases.  The  flow  indicated  by  the  rating  for  a  particular  stage  is  the  apparent  streamflow. 
Deviations  in  the  average  relationship  between  stage  and  streamflow,  which  are  interpreted  from 
individual  direct  measurements  of  them,  are  described  by  shifts  in  the  rating.  Together  with  hourly 
or  more  frequent  measurements  of  stage,  the  rating  and  its  shifts  are  used  to  compute  streamflow 
records  for  armual  publication,  following  a  comprehensive  analysis  of  streams  in  a  network  of 
gaging  stations. 

During  open-water  conditions,  uncertainty  associated  with  shifts  in  the  rating  usually  is  small, 
and  doesn't  much  affect  dam  operations.  Thus,  real-time  estimates  of  flow,  traditionally  deter¬ 
mined  directly  from  the  rating  and  in  some  cases  on  the  basis  of  a  shift  in  the  rating  defined  from  the 
most  recent  direct  measurement,  provide  sufficient  accuracy.  However,  during  ice-affected  periods, 
the  accuracy  of  real-time  streamflow  estimates  is  compromised  by  rapidly  varying  ice-backwater 
conditions  and  large  time- varying  shifts.  In  these  conditions,  the  ratio  of  the  true  to  the  apparent 
streamflow,  called  the  streamflow  ratio,  can  vary  from  1  during  open-water  conditions  to  a  value 
near  0.  These  shifts  create  enough  uncertainty  in  traditional  real-time  flow  estimates  that  dam  oper¬ 
ations  can  be  hindered. 

Previous  work 

Melcher  and  Walker  (1992)  evaluated  and  compared  17  methods  for  estimating  ice-affected 
streamflow  that  have  traditionally  been  or  could  potentially  be  used  in  the  nationwide  streamflow- 
gaging  station  network  maintained  by  USGS.  The  methods  were  divided  into  two  general  cate¬ 
gories,  subjective  and  analytical,  depending  on  whether  or  not  judgment  was  necessary  for  method 
application.  They  identified  two  subjective  methods  that  were  more  accurate  than  other  subjective 
methods  analyzed,  and  approximately  as  accurate  as  the  best  analytical  method.  Holtschlag  (1996) 
developed  a  dynamical-systems  approach  for  computing  ice-affected  streamflow.  This  approach 
ranked  higher  than  the  11  analytical  methods  investigated  by  Melcher  and  Walker  (1992)  on  the 


basis  of  accuracy  and  feasibility  criteria.  The  difference  equation  formulation  developed  by 
Holtschlag  (1996)  provided  the  basis  for  the  filter  presented  in  this  report. 

Purpose  and  scope 

This  report  describes  an  extended  Kalman  filter  that  can  be  used  to  project  the  expected  value 
and  the  associated  uncertainty  of  ice-affected  streamflow.  The  filter  was  applied  and  evaluated  on 
the  basis  of  historical  climatological  and  streamflow  data  available  near  the  St.  John  River  at  Dickey, 
Maine,  and  near  the  Platte  River  at  North  Bend,  Nebraska.  The  filter  was  restricted  to  forms  that 
could  be  automated  and  applied  by  use  of  data  that  have  been  routinely  compiled  and  that  could  be 
readily  obtained  in  real  time. 


SITES  OF  DATA  COLLECTION 

Two  sites  for  development  and  initial  application  of  the  filter  were  selected  to  represent  a  sample 
of  severe  ice-backwater  conditions  in  the  U.S.  In  addition,  the  selected  USGS  streamflow  gaging 
stations,  St.  John  River  at  Dickey,  Maine  (station  01010500),  and  Platte  River  at  North  Bend, 
Nebraska  (station  06796000),  have  a  long  historical  record  of  streamflow  data  and  are  located  close 
enough  to  climatological  data  collection  sites  to  support  filter  development  and  evaluation. 

St.  John  River  at  Dickey,  Maine,  gage 

The  gage  is  located  in  northern  Maine,  about  160  miles  north  of  Bangor,  and  about  9  miles  south¬ 
west  of  the  border  with  Canada  (Fig.  1).  It  is  situated  on  the  south  (right)  bank  of  the  river,  about  500 
ft  downstream  from  the  bridge  over 
the  St.  John  River  near  the  end  of 
State  Route  161.  The  drainage  area  at 
the  gage  is  2680  mi^. 

Daily  streamflow  data  have  been 
compiled  at  the  gaging  station  since 
September  1946.  Station  records 
through  water  year*  (WY)  1995  show 
a  maximum  instantaneous  flow  of 
91,700  ft^/s  on  29  April  1979,  a  mini¬ 
mum  daily  mean  flow  of  135  ft^/s  on 
15  September  1948,  and  an  average 
of  4771  ftVs.  The  stage-streamflow 
relation  is  usually  affected  by  ice 
from  early  December  through  mid- 
April.  Minimum  and  maximum  air 
temperature  data  are  recorded  at  the 
Fort  Kent  climatological  station, 
which  is  located  about  25  miles  east- 
northeast  of  the  gaging  station  (Fig. 

1).  Climatological  data  for  the  Fort 
Kent  station  are  published  by  the 
National  Oceanic  and  Atmospheric 
Administration  under  station  num¬ 
ber  2878.  ■  Figure  1.  Selected  streamfloio-gaging  and  climatological  sta¬ 

in  this  report,  17  years  of  daily  tions  near  the  St.  John  River  at  Dickey,  Maine. 

*  A  water  year  is  the  12  months  from  1  October  to  30  September  of  the  next  year.  For  example,  WY  1990  runs  from  1  October 
1989  thorough  30  September  1990. 
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Figure  2.  Distribution  of  streamflow  ratios  inferred  from  selected  ice-affected  periods 
at  St.  John  River  at  Dickey,  Maine,  gage  from  1971  through  1993. 


streamflow  and  climatolological  data  were  analyzed.  The  periods  included  WY 1971-74, 1978-79, 
and  1983-93;  intervening  periods  were  omitted  owing  to  intervals  of  missing  stage  or  climatololog¬ 
ical  data.  Within  the  selected  periods,  37.1%  (2303  days)  of  the  daily  values  were  from  ice-affected 
conditions,  as  shown  by  streamflow  ratios  of  less  than  1.  The  distribution  of  streamflow  ratios, 
inferred  by  dividing  published  values  by  corresponding  apparent  values,  is  right  skewed,  with  a 
mode  class  of  0.15  (Fig.  2).  For  the  selected  periods,  the  minimum  streamflow  ratio  was  0.014  and 
the  mean  of  ice-affected  streamflow  ratios  was  0.23.  The  corresponding  distribution  of  daily  mean 
air  temperature  values,  computed  as  the  average  of  the  daily  maximum  and  minimum  values,  on 
days  of  ice  effects  is  approximately  symmetrically  distributed  with  a  mean  of  -9.4°C  (Fig.  3). 


Average  Air  Temperature  (°C) 

Figure  3.  Distribution  of  average  daily  air  temperatures  at  Fort  Kent,  Maine,  during 
selected  ice-affected  periods  at  St.  John  River  at  Dickey,  Maine,  gage  from  1971  through 
1993. 
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Figure  4.  Selected  streamflow-gaging  and  climatic  sta¬ 
tions  near  the  Platte  River  at  North  Bend,  Nebraska. 


Figure  5.  Distribution  of  streamflow  ratios  inferred  from  selected  ice-affected  periods 
at  Platte  River  at  North  Bend,  Nebraska,  gage  from  1965  through  1994. 
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Figure  6.  Distribution  of  average  daily  air  temperatures  at  Fremont,  Nebraska, 
during  selected  ice-affected  periods  at  Platte  River  at  North  Bend,  Nebraska,  gage 
from  1965  through  1994. 

Platte  River  at  North  Bend,  Nebraska,  gage 

This  gage  is  located  in  east-central  Nebraska,  about  45  miles  northwest  of  Omaha.  It  (USGS 
gaging  station  06796000)  is  situated  on  the  north  (left)  bank  of  the  river,  about  80  ft  west  (upstream) 
of  State  Highway  79  (Fig.  4)  and  about  0.7  miles  south  of  North  Bend,  Nebraska.  The  drainage  area 
of  the  gage  is  70,400  mi^.  Daily  streamflow  data  have  been  compiled  at  the  gaging  station  since 
April  1949.  Station  records  through  WY 1995  show  a  maximum  instantaneous  flow  of  112,000  ft^/s 
on  29  March  1960,  a  minimum  daily  mean  flow  of  36  ft^/ s  on  29  July  1974,  and  an  average  of  4569 
ft^/s.  Air  temperature  data  were  obtained  from  a  weather  station  at  Fremont,  Nebraska,*  which  is 
located  about  15  miles  east  of  the  gaging  station  (Fig.  4). 

Streamflow  at  the  gaging  station  is  subject  to  a  complete  ice  cover  during  the  winter.  Ice  effects 
are  common  in  November  through  March;  there  are  severe  ice  jams  at  the  gage  in  most  years.  In  this 
report,  27  years  of  daily  streamflow  and  climatological  data  were  analyzed:  WY  1965-67, 1969-70, 
1972-90,  and  1992-94.  Intervening  periods  were  omitted  owing  to  missing  stage  or  climatological 
data. 

Within  the  period  analyzed,  28.9%  (2179  days)  of  daily  streamflow  values  were  affected  by  ice,  as 
indicated  by  a  streamflow  ratio  of  less  than  1.  The  distribution  of  streamflow  ratios  during  periods 
of  ice  effects  is  skewed  to  the  right,  with  a  mode  class  of  0.30  (Fig.  5).  For  the  periods  of  analysis,  the 
minimum  streamflow  ratio  was  0.045  and  the  mean  of  ice-affected  streamflow  ratios  was  0.43.  The 
corresponding  distribution  of  daily  mean  air  temperatures  on  days  of  ice  effects  has  a  mode  class  of 
0°C  and  a  mean  of  -4.1°C  (Fig.  6). 


APPLICATION  OF  EXTENDED  KALMAN  FILTERING 

Grewal  and  Andrews  (1993),  Bar-Shalom  and  Li  (1993),  Bozic  (1994),  Mendel  (1995),  and  Brown 
and  Hwang  (1997)  all  provide  detailed  information  on  the  mathematical  development  and  general 
application  of  Kalman  filtering.  A  Kalman  filter  estimates  the  state  of  a  dynamic  system,  given 
measurements  that  are  related  to  the  state.  The  extended  form  of  the  Kalman  filter  was  selected 

*  Personal  communication  with  Mat  Werner,  Climate  Resources  Specialist,  High  Plains  Climate  Center,  University  of 
Nebraska-Lincoln,  1996. 
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because  the  state  vector  was  formulated  to  include  both  a  signal  element  (the  streamflow  ratio)  and 
unknown  parameters.  In  addition  to  the  imknown  parameters  in  the  state  vector,  five  threshold 
parameters  were  estimated  externally  to  the  extended  Kalman  filter. 

A  discrete-time  formulation  was  used  for  consistency  with  the  availability  of  hydrological  and 
climatological  data  and  for  filter  simplicity.  In  this  report,  the  length  of  the  discrete  time  step  was 
1  day,  a  length  that  facilitated  analysis  of  extensive  historical  periods  of  record. 

A  discrete-time  extended  Kalman  filter  was  developed  to  accoimt  for  the  effects  of  ice  on  stream- 
flow.  The  filter  consists  of  two  models,  a  nonlinear  process  model  and  a  linear  measurement  model. 
The  general  form  of  the  nonlinear  process  model  is 

xik)  =  f[x(k-l),k-l]  +  w(k-l)  (1) 

where  x(fc)  =  state  vector.  In  this  report,  the  state  vector  is  partitioned  into  two  components.  The 
first  is  the  streamflow  ratio  and  is  called  the  state  signal  element.  The  second  is  the 
state  parameter  vector.  The  total  number  of  elements  in  the  state  vector  is  the  dimen¬ 
sion  of  the  state  space. 

/[x(fc-l),  fc-1]  =  nonlinear  function  of  the  state  at  the  previous  time  step  plus  other  information  on 
auxiliary  variables  available  at  time  k-1. 

w(/c-l)  =  value  from  a  random  sequence  representing  process  noise*  at  time  k-1.  The 
sequence  w  is  assumed  to  be  independent  and  normally  distributed,  with  a  mean  of 
zero  and  a  diagonal  covariance  structure  Q{k-1)  commonly  written  w~N[0,  Q(A:-1)]. 
In  this  application,  only  the  variance  of  Xi{k)  was  assumed  be  nonzero;  no  process 
noise  was  associated  with  the  state  parameters. 

The  time-varying  linear  measurement  model  is  of  the  form 

z{k)  =  H{k)  \{k)  +  v(k)  (2) 

where  z(k)  =  streamflow  at  time  k. 

H(k)  =  time-varying  measurement  sensitivity  matrix  that  is  represented  by  the  vector 
lh(k)  0  0],  where  h(k)  is  the  apparent  streamflow. 
v(k)  =  value  from  a  random  sequence  representing  measurement  noise  at  time  k.  The 
sequence  v  is  assumed  to  be  independent  and  normally  distributed,  with  a  mean  of 
zero  and  variance  of  R{k).  The  subset  of  days  indexed  by  k  on  which  direct  measure¬ 
ments  of  streamflow  were  obtained  is  denoted  as  V.  For  developing  a  projection,  the 
variance  R{k')  was  assumed  to  be  proportional  to  the  published  streamflow  on  days 
of  direct  measurement.  In  open-water  conditions,  the  standard  error  was  assumed  to 
be  2.5%  of  the  published  values;  during  ice-affected  conditions,  the  standard  error 
was  assumed  to  be  8.0%  of  the  published  streamflow.  On  days  without  direct  meas¬ 
urement,  published  flow  data  were  not  used  to  update  the  estimate  of  the 
streamflow  ratio. 

Formulation 

Based  upon  analysis  of  historical  streamflow  characterizations  developed  by  traditional  meth¬ 
ods  for  estimating  ice-affected  flow,  the  dynamics  of  ice  effects  were  classified  into  three  modes. 
Mode  1  dynamics  are  generated  by  the  sudden  formation-ablation  of  ice,  as  indicated  by  abrupt 
changes  in  apparent  streamflow.  Mode  2  dynamics  are  caused  by  stable  ice  conditions  and  are 
approximated  by  a  first-order  difference  equation  relating  the  streamflow  ratio  to  air  temperature 


*  In  the  analysis  of  dynamic  systems,  noise  is  random  inputs  that  cannot  be  directly  measured  or  controlled. 


6 


values.  Mode  3  dynamics,  finally,  come  from  the  eventual  elimination  of  ice  effects  at  higher  air 
temperatures.  Each  mode  corresponds  to  a  xmique  process  model;  how^ever,  the  measurement 
model  was  consistent  for  each  mode. 

The  process  model  for  mode  1  dynamics  is  an  algebraic  equation  of  the  form 

=  (3) 

Mode  1  dynamics  were  applied  for  either  of  two  conditions: 

•If  the  1  day  change  in  apparent  streamflow  increased  by  more  than  q_dl  percent  and  the  average 
air  temperature  was  less  then  t_lo. 

•If  the  one  day  change  in  apparent  streamflow  decreased  by  more  than  q_dl  percent  when  the 
streamflow  ratio  was  less  than  1. 

Because  of  the  lack  of  derivative  information  in  eq  3,  q_dl  and  t_lo  were  included  among  five  thresh¬ 
old  parameters  estimated  outside  the  extended  Kalman  filter.  Thus,  for  mode  1  dynamics,  the  state 
vector  only  included  the  streamflow  ratio  and  the  corresponding  dimension  of  the  state  space  was  1. 

The  process  model  for  mode  2  dynamics  was  a  first-order  difference  equation,  driven  by  air 
temperatures  u,  that  was  similar  to  one  developed  and  implemented  in  a  previous  investigation 
(Holtschlag  1996).  The  process  model  for  mode  2  dynamics  is  of  the  form 

^1  (fc)  =  ^2  (fc  - 1)  +  ^3  [3^1  (fc  - 1)  -  X2  (k  - 1)]+  Xi(k  -  l)[u{k  - 1)  -  Xs  (*:  - 1)]  (4) 

which  says  that  at  times  of  ice  effects  and  constant  air  temperature  x^,  the  streamflow  ratio  (xj)  is  in 
equilibrium  about  a  nominal  value  X2-  Changes  from  the  equilibrium  value  are  described  by  a  dif¬ 
ference  equation  that  includes  an  autoregressive  component  with  parameter  x^  and  a  forcing  fimc- 
tion  term  driven  by  daily  air  temperature  values.  Air  temperatures  that  vary  from  a  nominal  value 
of  X5  change  the  streamflow  ratio  from  its  nominal  value  of  3:2  at  a  rate  of  X4.  This  form  of  a  difference 
equation  is  nonlinear  in  parameters  because  both  rate  parameters  (X3  and  X4)  and  offset  parameters 
{x2  and  Xs)  are  estimated  simultaneously. 

Prior  information  on  the  distribution  of  streamflow  ratios  during  periods  of  ice  effects  was  used 
in  hope  of  facilitating  the  solution  of  this  inherently  difficult  estimation  problem.  Air  temperature 
values  u  indicate  an  exponentially  weighted  average  of  temperatures  from  the  3  previous  days.  An 
exponential  weighting  factor  t_wt,  used  in  computing  u,  was  included  among  threshold  para¬ 
meters,  because  it  did  not  occur  explicitly  in  eq  4.  The  dimension  of  the  state  space  for  mode  2 
dynamics  is  5. 

The  process  model  for  mode  3  dynamics  is  an  algebraic  equation  of  the  form 

xi(fc)  ^  - 1) + y  ~  ~  [1  -  -1)].  (5) 

t_ou  -  t_hi 

Mode  3  dynamics  are  applied  when  the  exponentially  weighted  air  temperature  value  exceeds  t_hi. 
Then,  the  streamflow  ratio  increases  from  its  value  at  time  A:  -  1  to  a  value  of  1  when  the  exponen¬ 
tially  weighted  air  temperatures  equal  t_ou.  Because  tjii  and  t_ou  occur  explicitly  in  the  process 
model,  they  could  be  included  as  parameters  in  the  state  vector.  However,  for  simplicity  in  this 
analysis,  they  were  included  among  threshold  parameters  estimated  outside  of  the  extended  Kal¬ 
man  filter.  As  in  the  case  of  mode  1  dynamics,  the  dimension  of  the  state  space  for  mode  3  dynamics 
was  treated  as  1.  Estimates  of  streamflow  ratio  for  all  dynamic  modes  were  constrained  to  the  inter¬ 
val  between  the  minimum  ratio  of  published  to  apparent  streamflow  determined  from  historical 
record  and  1. 

The  process  models  for  the  three  dynamic  modes  were  developed  so  that  the  first  (only)  element 
in  the  state  vector  was  the  signal  element  (the  estimated  streamflow  ratio).  For  this  application. 
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changes  in  mode  only  affected  the  signal  element.  In  contrast,  the  state  parameter  vector  was  treat¬ 
ed  as  if  mode  2  dynamics  were  always  taking  place.  Although  this  convention  can  create  imcer- 
tainty  in  the  state  parameter  vector,  the  effect  is  lessened  because  updates  affecting  the  parameter 
vector  are  only  computed  for  days  of  direct  streamflow  measurement.  Direct  measurements  are  not 
possible  during  mode  1  or  3  dynamics  because  of  unsafe  measuring  conditions. 

Implementation 

The  extended  Kalman  filter  was  implemented  by  recursively  computing  daily  updates  to  the 
state  vector  x  and  the  state  error  covariance  matrix  P,  given  an  initial  estimate  of  the  state  vector  xq 
and  Pq.  The  update  has  two  steps:  a  temporal  update  and  an  observational  update.  The  temporal 
update  is  computed  at  each  time  step;  the  observational  update  is  computed  only  on  days  of  direct 
measurement.  The  magnitude  of  P  increases  with  temporal  updates  and  decreases  with  observa¬ 
tional  updates.  Bootstrap  estimates  of  xq  and  Pq  were  computed  by  iteratively  replacing  x|/'*’^^  with 
until  where  Xq  and  Xf  are  the  initial  and  final  values  for  the  state  vector  for  the  j  or 

;  +  1  iteration  of  the  extended  Kalman  filter,  respectively. 

Temporal  updates 

The  temporal  update  represents  the  best  linear  estimate  of  the  state  at  time  k  without  consider¬ 
ation  of  available  observations  at  time  k,  but  only  information  available  at  k-1.  The  notational  form 
used  for  the  general  state  equation  is 

x(-)(fc)=/{xW[(fc-l),fc-l]}.  (6) 

For  mode  1  dynamics,  the  signal  element  is  replaced  by 

x[-Hk)  =  ^^^^.x,{k-l).  (7) 

In  the  case  of  mode  2  dynamics,  the  temporal  update  is  written 


'x[-\k) 

'x^^Hk  - 1)  -t  - 1)  -  - 1)]  +  -  l)[u{k  - 1)  -  - 1)]' 

x[^\k-l) 

x^i\k) 

= 

1 

xi^Hk-1) 

x^iHk\ 

x^^\k-V) 

Finally,  in  the  case  of  mode  3  dynamics,  the  signal  element  is  computed  as 


x\-\k)  =  x{'^Hk-l)  + 


u{k-T)-t_hi 

t_ou-t_hi 


(9) 


for  specified  threshold  parameters  tjii  and  t_ou. 

A  temporal  update  of  streamflow  is  computed  by  multiplying  the  apparent  streamflow  by  the 
estimated  streamflow  ratio,  or,  for  consistency  with  the  extended  Kalman  filter  notation,  by  multi¬ 
plying  the  time-varying  measurement  sensitivity  matrix  H{k)  times  the  temporal  update  of  the  state 
vector  as 

zi-\k)  =  H{k)xi-'i{k)  =  [h{k)0  0]x(-'>{k).  (10) 

A  temporal  update  of  the  error  covariance  matrix  is  computed  as 
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(11) 


p(-)  (fc)  =  (A:  -  (fc  - 1)  (fc  - 1)  +  Q{k  - 1). 

A  first  order  approximation  of  the  state  transition  matrix  is  given  by 

f'dc-D-l/Ij:,*:-!!,.,, 

The  estimate  of  the  transition  matrix  used  in  this  analysis  is 

0  10  0  0 

Oh](A;-l)  =00  1  00 

0  0  0  1  0 

0  0  0  0  1 

The  value  of  Q  was  determined  such  that 

Prob[|(-)(fc)  ±  Za=o.i  Pl:^ik)  >  z{k)]  ~  0.9  (14) 

where  Z„=o,i  is  the  standardized  normal  deviate  corresponding  to  a  90%  probability  and  equaling 
1.64. 

Observational  updates 

Observational  updates  were  computed  for  days  of  direct  flow  measurement.  First,  the  Kalman 


gain  matrix  was  computed  as 

Kik')  =  Pi-Hk')HT{k')[H{k')PO)(k')HT{k')  +  R{k')Y\  (15) 

Then  the  covariance  matrix  was  updated  as 

P^^\k')  =  [/  -  K{k')Hik')]P^-\k').  (16) 

The  state  vector  update  was  computed  as 

x^*\k')  =  x^-\k')  +  K{k')[z{k')  -  z^-\k')].  (17) 

And  finally  the  observational  update  was  computed  as 

zW(fc)  =  H{k)x^^\k)  =  [iz(fc)00]x^+)(fc).  (18) 


On  days  without  direct  flow  measurement,  observational  updates  were  set  equal  to  the  temporal 
updates  computed  for  that  time  step.  Thus,  projected  values  computed  by  use  of  the  extended  Kal¬ 
man  filter  were  equal  to  the  observational  updates  on  days  of  direct  streamflow  measurement  and 
were  equal  to  the  temporal  updates  otherwise.  No  adjustment  was  included  for  uncertainty  in  the 
apparent  streamflow  values. 

Computation 

Computational  procedures  were  developed  and  tested  in  the  MATLAB  programming  environ¬ 
ment  (The  Math  Works,  Inc.,  Natick,  Massachusetts)  on  the  basis  of  algorithms  developed  by  Grewal 
and  Andrews  (1993).  To  provide  numerical  stability,  the  temporal  update  of  the  error  covariance 
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matrix  was  computed  by  use  of  the  Thornton  UD  factorization  algorithm  (Grewal  and  Andrews 
1993,  p.  255)  rather  than  by  directly  using  eq  6-13.  Similarly,  the  Bierman  observational  update  algo¬ 
rithm  (Grewal  and  Andrews  1993,  p.  245)  with  corrections  provided  by  D.  Flynn*  was  used  rather 
than  directly  using  eq  15-17.  The  MATLAB  code  was  converted  to  code  in  the  C-I-+  programming 
language  (Weiss  1996)  for  improved  portability  (a  program  listing  is  available  from  the  first  author, 
DJH). 

PROJECTING  ICE-AFFECTED  STREAMFLOW 

Filter  estimates  generally  are  at  time  k,  based  on  measurements  up  to  and  including  time  k.  In 
contrast,  forecast  estimates  are  made  at  time  k,  based  on  data  up  to,  but  not  including,  time  k.  In  this 
report,  extended  Kalman  filter  estimates  of  streamflow  are  projections,  forecasts  determined  on  the 
basis  of  the  temporal  updates.  However,  on  days  of  direct  measurement,  more  accurate  filter  esti¬ 
mates  are  computed  by  use  of  observational  updates. 

St.  John  River  at  Dickey,  Maine 

The  extended  Kalman  filter  was  initialized  to  St.  John  River  data  by  manually  adjusting  prelimi¬ 
nary  estimates  of  threshold  parameter  values.  This  minimized  the  sum  of  squared  errors  in  the 
extrapolated  streamflow  ratio  (k')  -  Xi(k') ,  where  k'  indicates  days  of  direct  streamflow  measure¬ 
ments.  Once  satisfactory  estimates  of  the  threshold  parameters  were  obtained,  they  were  fixed  (Table 
1).  Then  the  filter  was  run  repetitively,  using  the  state  vector  and  error  covariance  matrix  computed 
on  the  last  iteration  of  the  previous  run  to  initialize  the  subsequent  nm.  The  filter  was  run  repeatedly 
imtil  elements  in  the  state  parameter  vector  were  essentially  constant.  In  this  process,  initial  esti¬ 
mates  for  the  state  error  covariance  matrix  converged  from  an  initially  specified  diagonal  matrix 
with  nonzero  components  of  [0.5  0.5  0.01  10]  to  the  standard  errors  shown  in  Table  1. 

Final  estimates  for  the  state  parameters  indicate  that  mode  1  dynamics  are  highly  autoregressive, 
as  shown  by  the  parameter  X3  =  0.981,  about  a  streamflow  ratio  offset  of  X2  =  0.544.  Streamflow  ratios 
increase  at  a  rate  x^  -  0.000855°C“l  from  the  temperature  offset  X5  =  -3.19°C.  Although  the  value  for  X2 


Table  1.  Threshold  and  filter  parameters  for  the  extended  Kalman  filter  for  the  St.  John 
River  at  Dickey,  Maine. 


Parameter 

symbol 

Parameter  description 

Estimate 

Sensitivity 

Estimated 

standard 

error 

tjii 

High  temperature  at  which  ice  breakup  begins 
(in  degrees  Celsius). 

4.50 

-0.0824 

_ 

tjo 

Low  temperature  at  which  sudden  increases  in  apparent 
streamflow  indicates  ice  accumulation  (in  degrees  Celsius). 

-2.25 

0.00868 

— 

t_ou 

High  temperature  at  which  ice  breakup  is  complete 
(in  degrees  Celsius). 

10.0 

0.0427 

— 

t_wt 

Exponential  weighting  factor  for  daily  temperatures. 

0.95 

0.0424 

— 

qjl 

Threshold  at  which  changes  in  apparent  daily  streamflows 
are  considered  large. 

0.35 

-16.1 

_ 

Xl 

Offset  streamflow  ratio. 

0.544 

— 

0.0937 

Autoregressive  parameter  for  streamflow  ratio. 

0.981 

— 

0.000153 

3:4 

Parameter  relating  air  temperature  to  changes  in 
streamflow  ratios  (in  degrees  Celsius"'). 

0.000855 

_ 

0.00000459 

X5 

Offset  temperature  (in  degrees  Celsius). 

-3.19 

— 

2.08 

*  Personal  communication,  California  State  University  at  Fullerton,  1996. 
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is  higher  than  the  mode  class  of  the  distribution  of  empirical  streamflow  ratios  (Fig.  2),  it  is  physi¬ 
cally  realizable.  Similarly,  X5  is  consistent  with  the  distribution  of  air  temperatures  during  periods 
of  ice  effects  and  the  physical  conceptualization  of  the  process  (Fig.  3).  However,  analysis  of  the 
state  error  covariance  matrix  points  to  a  large  positive  correlation  (>0.99)  between  X2  and  X5  and  a 
large  negative  correlation  (<  -0.93)  between  3:3  and  X4.  Thus,  although  the  form  of  the  difference 
equation  used  to  describe  mode  1  dynamics  resulted  in  parameters  with  physically  realizable  val¬ 
ues,  the  large  correlations  indicate  ambiguity  concerning  their  true  values.  Because  of  the  high 
correlations  in  the  state  error  covariance  matrix,  there  is  a  potential  for  reducing  the  dimension  of 
the  state  parameter  vector  without  loss  of  filter  accuracy. 

Sensitivities  for  threshold  parameters  (Table  1)  were  estimated  as  the  change  in  the  sum  of 
squared  errors  in  the  streamflow  ratio  divided  by  the  change  in  the  corresponding  parameter  near 
the  selected  values.  Results  of  simulations  show  that  filter  computations  were  most  sensitive  to 
changes  in  the  q_dl  parameter  and  least  sensitive  to  the  t_lo  parameter.  Formal  optimization  of  the 
threshold  parameters  could  lead  to  further  improvement  in  filter  performance. 

One  measure  of  the  projection  accuracy  of  the  filter  is  the  relation  between  temporal  updates  of 
streamflow  z[~\k')  and  published  daily  values  z(fc'),  where  fc'  indexes  days  of  direct  measurement 
(Fig.  7).  Although  this  value  is  modified  by  observational  updates  to  more  precisely  project  stream- 
flow  on  days  of  direct  measurement,  z[~\k')  provides  a  conservative  indication  of  the  filter  accu¬ 
racy.  The  temporal  update  is  a  conservative  sign  of  accuracy  because  the  variance  of  the  projection 
increases  monotonically  with  time  from  the  previous  measurement,  and  the  length  of  the  projection 
is  at  maximum  just  before  the  observational  update.  Results  for  the  St.  John  River  show  that  the 
correlation  between  log-transformed  values  of  z^\k')  and  ziV),  based  on  40  days  of  ice-affected 
measurements,  is  0.777  and  is  0.998  based  on  63  days  of  open-water  measurements.  Measurements 
at  the  St.  John  River  used  in  this  analysis  averaged  8.6  weeks  apart. 

Another  measure  of  filter  accuracy  is  the  relationship  between  published  and  projected  stream- 
flows  during  periods  of  ice  effects.  This  relation  is  linear  in  the  logarithm  of  transformed  values 
(Fig.  8).  Uncertainty  occurs  in  both  the  published  and  projected  values.  Published  values  during  ice 
effects  are  subjectively  rated  "fair"  or  "poor."  A  rating  of  "fair"  implies  that  about  95%  of  the  daily 
values  are  within  15%  of  the  true  value;  a  "poor"  rating  means  that  daily  streamflow  values  have 


Figure  7.  Relationship  between  published  streamflow  and  temporal  updates 
of  streamflow  on  selected  days  of  direct  measurement  at  St.  John  River  at 
Dickey,  Maine,  gage  from  1971  through  1993. 
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less  than  "fair"  accuracy  (Novak  1985,  p.  65). 
Analysis  of  the  distribution  of  discrepancies  (Fig. 
9)  between  published  and  projected  values  during 
periods  of  ice  effects,  computed  as 

^  log[z(+)(fc)]-log[2(fc)] 
log[z(lc)] 

indicates  that  the  absolute  value  of  elements  in  the 
e  sequence  are  less  than  8%,  87.2%  of  the  time  and 
are  less  than  15%,  96.6%  of  the  time. 

Projections  of  streamflow  are  shown  with  other 
streamflow  and  climatological  data  for  the  selected 
periods  in  Appendix  A.  Upper  and  lower  projec¬ 
tions  were  computed  by  adjusting  the  variance  of 
Q  to  0.0035  so  that  the  interval  formed  between  the 
upper  and  lower  projections  about  the  temporal 
updates  at  k’  contained  the  published  values  90% 
of  the  time. 


Figure  8.  Relationship  between  published  and 
projected  streamflow  during  selected  ice-affected 
periods  at  St.  John  River  at  Dickey,  Maine,  gage 
from  1971  through  1993. 


Platte  River  at  North  Bend,  Nebraska 

The  protocol  developed  for  the  St.  John  River  was  followed  when  the  extended  Kalman  filter 
was  applied  to  Platte  River  data.  This  was  done  by  manually  adjusting  preliminary  estimates  of 
threshold  parameters  to  minimize  the  sum  of  squared  errors  in  the  temporal  updates  of  the  stream- 
flow  ratio  for  days  of  direct  measurement.  Again,  once  apparently  satisfactory  estimates  of  the 
threshold  parameters  were  obtained,  they  were  fixed  (Table  2)  and  the  filter  was  nm  until  conver¬ 
gence.  In  this  process,  initial  estimates  of  the  state  error  covariance  matrix  converged  from  an  initial 
diagonal  matrix  with  nonzero  elements  of  [0.5  0.5  0.01  10]  to  a  covariance  matrix  with  standard 
errors  shown  in  Table  2. 


Discrepancy  Between  Logarithms  of  Published  and  Projected  Streamflow 


Figure  9.  Distribution  of  discrepancies  between  logarithms  of  published  and  pro¬ 
jected  streamflow  during  selected  ice-affected  periods  at  St.  John  River  at  Dickey, 
Maine,  gage  from  1971  through  1993  (logarithms  are  in  base  10). 
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Table  2.  Threshold  and  filter  parameters  for  the  extended  Kalman  filter  for  the  Platte  River  at 
North  Bend,  Nebraska. 


Parameter 

symbol 

Parameter  description 

Estimate 

Sensitivity 

Estimated 

standard 

error 

t_hi 

High  temperature  at  which  ice  breakup  begins  (in  degrees 
Celsius). 

4.50 

-0.331 

— 

tJo 

Low  temperature  at  which  sudden  increases  in  apparent 
streamflow  indicates  ice  accumulation  (in  degrees  Celsius). 

-2.25 

-0.141 

— 

t_ou 

High  temperature  at  which  ice  breakup  is  complete  (in 
degrees  Celsius). 

10.0 

-0.206 

— 

t_wt 

Exponential  weighting  factor  for  daily  temperatures. 

0.95 

0.430 

— 

qjl 

Threshold  at  which  changes  in  apparent  daily  streamflows 
are  considered  large. 

0.30 

10.5 

— 

Xj 

Offset  streamflow  ratio. 

-0.0681 

— 

0.0132 

xz 

Autoregressive  parameter  for  streamflow  ratio. 

0.990 

— 

0.000186 

Xi 

Parameter  relating  air  temperature  to  changes  in  streamflow 
ratios  (in  degrees  Celsius-'). 

0.000939 

— 

0.00000372 

X5 

Offset  temperature  (in  degrees  Celsius). 

-9.37 

— 

0.169 

Results  from  Platte  River  data  also  point  out  that  mode  1  dynamics  are  highly  autoregressive,  as 
indicated  by  the  parameter  jc3=0.990.  Streamflow  ratios  increase  at  a  rate  X4  =  0.000939°C“^  about  a 
temperature  offset  X5  =  -9.37°C,  a  lower  temperature  than  that  estimated  for  the  St.  John  River. 
Unfortunately,  the  estimated  streamflow  ratio  offset  of  X2  =  -0.068  is  not  physically  realizable.  Anal¬ 
ysis  of  the  state  error  covariance  matrix  shows  a  maximum  positive  correlation  of  0.81  between  X3 
and  X5  and  a  maximum  negative  correlation  of  -0.74  between  X3  and  X4.  The  magnitudes  of  these 
correlations  are  not  thought  to  be  sufficient  to  significantly  degrade  parameter  estimates.  However, 
given  the  small  magnitude  of  the  estimated  X2  value,  in  future  applications  it  may  be  possible  to 
eliminate  (set  to  zero)  the  streamflow  ratio  offset  from  the  difference  equation  for  mode  1  dynamics. 
Such  an  elimination  would  reduce  the  dimension  of  the  state  space,  which  would  also  likely  reduce 
parameter  ambiguity  caused  by  high  correlations  in  the  state  error  covariance  matrix.  Values  for  the 
threshold  parameters  that  are  less  than  optimal  also  possibly  explain  the  discrepancy  between  the 
estimated  value  of  X2  and  the  conceptualized  value. 

Sensitivities  for  threshold  parameters  (Table  2)  were  estimated  as  the  change  in  the  sum  of 
squared  errors  in  the  streamflow  ratio  estimate  divided  by  the  change  in  the  corresponding  param¬ 
eter  near  the  selected  values.  The  results  of  simulations  indicate  that  projections  are  most  sensitive 
to  changes  in  the  q_dl  parameter  and  least  sensitive  to  the  tjo  parameter.  Again,  formal  optimiza¬ 
tion  of  the  threshold  parameters  could  lead  to  further  improvement  in  filter  performance. 

The  temporal  updates  of  streamflow  on  days  of  direct  measurements  compare  closely  with  pub¬ 
lished  daily  mean  values  (Fig.  10).  Results  show  that  the  correlation  between  log-transformed  val¬ 
ues  of  z[~\k')  and  z(fc'),  based  on  87  days  of  ice-affected  measurements,  is  0.864  and  is  0.997  based 
on  345  days  of  open-water  measurements.  Measurements  at  the  Platte  River  used  in  this  analysis 
averaged  3.2  weeks  apart. 

The  relationship  between  published  and  projected  streamflow  values  at  the  Platte  River  during 
periods  of  ice  effects  is  linear  and  unbiased  in  the  logarithms  of  streamflow  (Fig.  11).  The  distribu¬ 
tion  of  discrepancies  between  published  and  projected  values  (Fig.  12)  during  periods  of  ice  effects 
were  analyzed  by  use  of  eq  18;  the  absolute  value  of  elements  in  the  e  sequence  are  less  than  8%, 
90.7%  of  the  time,  and  are  less  than  15%,  97.7%  of  the  time. 

Projected  streamflows  are  shown  with  other  flow  and  climatological  data  for  the  selected  peri¬ 
ods  in  Appendix  A.  Upper  and  lower  projections  were  computed  by  adjusting  the  variance  of  Q  to 
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Figure  10.  Relationship  between  published  streamflow  and  temporal  updates  of 
streamflow  on  selected  days  of  direct  measurement  at  Platte  River  at  North  Bend, 
Nebraska,  gage  from  1965  through  1994. 


Figure  11.  Relationship  between  published  and  pro¬ 
jected  streamflow  during  selected  ice-affected  peri¬ 
ods  at  Platte  River  at  North  Bend,  Nebraska,  gage 
from  1965  through  1994. 

0.0039  so  that  the  interval  formed  by  the  upper  and  lower  projection  about  the  temporal  update  at 
k'  included  the  published  values  90%  of  the  time. 


SUMMARY  AND  CONCLUSIONS 

Ice-affected  streamflow  describes  the  effects  of  charmel  ice  on  the  relationship  between  stream 
stage  (water-surface  elevation)  and  flow.  Because  an  ice  cover  physically  blocks  the  channel  and 
increases  flow  resistance,  these  effects  result  in  lower  flows  than  would  be  expected  from  corre- 
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Discrepancy  Between  Logarithms  of  Published  and  Projected  Streamflow 


Figure  12.  Distribution  of  discrepancies  between  logarithms  of  published  and  pro¬ 
jected  streamflow  during  selected  ice-affected  periods  at  Platte  River  at  North 
Bend,  Nebraska,  gage  from  1965  through  1994. 

sponding  stages  during  open-water  conditions.  In  addition,  ice  effects  create  uncertainty  in  real¬ 
time  streamflow  estimates  that  are  needed  to  help  control  floods  and  facilitate  navigation. 

Three  dynamic  modes  of  ice  effects  were  identified  on  the  basis  of  historical  interpretations  of 
ice-affected  streamflow  records.  In  the  first  mode  come  rapid  changes  in  ice  effects,  associated  with 
ice  formation  and  ablation,  indicated  by  abrupt  changes  in  apparent  streamflow.  In  the  second 
mode,  existing  ice  effects  change  with  changes  in  air  temperature.  In  the  third  mode,  ice  effects  are 
eliminated  with  higher  air  temperatures.  Equations  for  these  dynamic  modes  were  developed 
within  a  discrete-time  extended  Kalman  filter  to  project  daily  values  and  xmcertainties  of  ice- 
affected  streamflow. 

The  filter  consists  of  two  models:  a  nonlinear  process  model  and  a  time-varying  linear  measure¬ 
ment  model.  For  the  dominant  mode  2  dynamics,  the  process  model  computes  the  ratio  of  actual  to 
apparent  flow  (streamflow  ratio)  by  use  of  a  first-order  difference  equation  driven  by  daily  air- 
temperature  values  and  four  filter  parameters.  During  periods  of  mode  1  or  mode  3  dynamics,  the 
difference  equation  is  replaced  by  an  algebraic  expression  to  estimate  the  streamflow  ratio  by  use  of 
five  threshold  parameters. 

The  measurement  model  projects  a  daily  mean  streamflow  on  the  basis  of  the  estimated  stream- 
flow  ratio  and  the  apparent  streamflow.  On  days  of  direct  measurement,  the  projection  is  computed 
and  the  covariance  matrix  is  updated  to  account  for  this  information.  The  utility  of  the  filter  was 
evaluated  by  applying  it  to  historical  data  from  two  long-term  stream-gaging  stations. 

The  filters  developed  for  stream-gaging  stations  on  the  St.  John  River  at  Dickey,  Maine,  and  the 
Platte  River  at  North  Bend,  Nebraska,  were  stable,  and  parameters  converged  for  both  stations, 
allowing  projections  of  ice-affected  streamflows.  Results  for  the  gaging  station  at  the  St.  John  River 
indicate  that,  during  periods  of  ice  effects,  logarithms  of  projected  streamflow  values  were  within 
8%  of  the  logarithms  of  published  values  87.2%  of  the  time  and  within  15%,  96.6%  of  the  time. 
Results  for  the  gaging  station  at  the  Platte  River  indicate  that  logarithms  of  projected  streamflow 
values  were  within  8%  of  the  logarithms  of  published  daily  values  90.7%  of  the  time  and  within 
15%,  97.7%  of  the  time.  The  correlations  between  temporal  updates  and  published  values  of  stream- 
flow  on  days  of  direct  measurement  are  0.777  and  0.998,  for  data  from  the  St.  John  River,  and  0.864 
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and  0.997  for  data  from  the  Platte  River.  Analysis  of  the  state  error  covariance  matrix  for  both  rivers 
indicates  that  some  reduction  in  the  number  of  parameters  associated  with  the  difference  equation 
formulated  for  mode  1  dynamics  is  possible. 

The  extended  Kalman  filter  developed  in  this  report  provides  a  basis  for  projecting  ice-affected 
streamflow  at  other  gaging  stations  by  adjusting  filter  parameters  to  site-specific  conditions.  The 
filter  can  project  daily  mean  flow  during  periods  of  ice  effects  by  use  of  real-time  climatological  and 
hydrological  data. 
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